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Abstract 

The aim of this paper is to give an account of some applications 
of pseudodifferential calculus for solving linear wave equations in the 
limit of high frequency /short wavelength waves. More specifically, on 
using as a benchmark the case of electromagnetic waves propagating 
in a cold isotropic slowly space- and time-varying plasma, it is shown 
that, in general, linear plasma waves are governed by pseudodifferential 
operators. Thereafter, the asymptotic techniques for solving the cor- 
responding pseudodifferential wave equations are presented with em- 
phasis on the paraxial propagation of Gaussian wave trains in a cold 
isotropic plasma. Finally, it is addressed the unicity of the disper- 
sion tensor in terms of which the considered asymptotic solutions are 
determined. 



1 Pseudodifferential wave equations: basic defini- 
tions and symbol calculus 

Let us consider a generic linear wave equation in the (abstract) form 

(D.^)(x) = f(x), (1) 

where D is a matrix of (linear) operators acting on a multi-component wave- 
field i/y(x) and f(x) an external drive (sources). Both if} and f are defined on 
the iV-dimensional real space, with x = (x 1 , . . . ,x N ) Cartesian coordinates, 
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whereas D is assumed to be a pseudodifferential (^DO) operator ^H], that 
is, its action on the wavefield xj) amounts to the Fourier integral operator 

(D • V) (x) = J e ik < x ~^d(k, x,x')- ij}(x')d N x'd N k, (2) 

where d(k, matrix of smooth functions, referred to as symbols, that 

depend on the wavevector k = (k\, . . . , fejy) as well as on the outcoming (x) 
and incoming (x 1 ) variables. Moreover, symbols are characterized by the 
estimate [Q 

9%dgd(k,z) = 0(5- m+ ^), z = (x,x') (3) 

for 6 = \kL\~ l — > (semiclassical limit), where the scalelength L character- 
izes the nonuniformity of the operator, being defined by d%d ~ d/L^I; in 
equation (j3J), a = (a%, . . . , aj^) and = . . . , 02n) are multi-indices of 
dimension N and 2N, respectively, with, e.g., <9j? = <9' Q /dk^ 1 ■ ■ ■ dk^ N and 
\a\ = X^i a «' whereas m is a real number specifying the order of both the 
symbol and the corresponding operator. The condition (J3J) is automatically 
satisfied if d is a polynomial in k with smooth x-dependent (matrix-valued) 
coefficients, namely, d(k,x,x') = d(k,x) = ^| a |< m d a (x)k a with m a pos- 
itive integer, the corresponding ^DO amounting to the partial differential 
operator (PDO) D = Yl\ a \<m(~ i)^d a (x)d" . It is worth noting that in this 
case the response (D • tp) (x) of the operator depends on the derivatives of 
the wavefield of order < m only, i.e., D is a local operator; in order to ac- 
count for nonlocal effects in the response, polynomials (of integer order m) 
have been generalized to symbols (of real order m). 

As a consequence of the estimate (jHJ), symbols can be exploited in order 
to deal with the corresponding operators (symbol calculus). With this aim 
it is convenient to perform, in equation (J2J), the coordinate transformation 
(x, x') — > (r = qx+px' , s = x—x') with g+p = 1, redefine the symbol accord- 
ing to d(k, r, s) = d(k, r+ps, r—qs) and expand with respect to s. As a result 
one gets the sequence of symbols d ( f' p \k,x) = J2\ a \=£ 12T 9^dfd(k,x,0) of 
decreasing order (= m — I, cf., equation ©) for £ = 0, 1, . . ., which can be 
resummed to give the new symbol 

+00 

d^,x)~^d^ p) (fc,x), (4) 

1=0 

where the "~" means that, for every integer n, 

order m — n. The new symbol @ is referred to as the reduced symbol of 



the operator since it depends on a single x-variable rather than on both 
x and x'. In the semiclassical limit 5 = \kL\~ 1 — > 0, the series (J3J is 
asymptotically convergent and for practical applications only the first two 
terms are significant. The advantage of the reduced symbol is that one can 
replace D by any ^DO E)( 9 ' p ) corresponding to d^ q ' p \k,qx + px'), which 
depends on (x,x') through the convex combination qx + px' only. More 
specifically, the remainder R = D — D^'P) is a ^DO of order m = -co, i.e., 
the corresponding symbol d — d^' p ^ satisfies the estimate (jHJ) for any order 
m, thus R can be neglected in the semiclassical limit 5 = \kL\~ 1 — ► 0, being 
formally of order <5 + °°. 

In particular, two specific forms of the reduced symbol play a key role 
in physical applications: 

• the te/it-reduced symbol d^ = d^ 1,0 \ 

• the Weyl symbol = d^a>a); 

the former is the simplest generalization of the symbol of a partial differ- 
ential operator (cf., section 2), whereas the latter properly accounts for the 
nonuniformity of the medium [2] (cf., section 3). 

For any asymptotic form it is possible to represent the composition of 
operators through a suitable product-rule for symbols on noting that, given 
two scalar ^DO A and B with symbols a and b, repectively, their composi- 
tion C = A B corresponds to the symbol c(k,x,x') = a^ 1 ' ' (k,x)b^ 0,1 > (k,x') 
where the first (second) factor has been represented by means of the left 
(right, i.e., q = 0, p = 1) symbol pQ; on applying the reduction procedure to 
c(k, x, x'), cf., equation (j3J, and converting the form of the symbols a*- 1 ' ) and 
5(0,1) ( a g a j n by means of (jlj with d replaced by a^'°\k, x) and &( 0,1 ) (&;, x') 
as appropriate), one gets, to lowest significant orders, 

(D (D,(L) -da (i) 



dki dx l 

i 



m b <w) + ±{ a <w) }b <.w)} + ... , (5b ) 



for the left and the Weyl symbol, respectively, with {a, b} = Jf- — Jf 
being the Poisson brackets in the x-k phase space. The product-rules © can 
be readily generalized to matrix-valued symbols (on replacing the product 
of functions by the product of matrices) and constitute a key tool in sym- 
bol calculus as they allow to perform calculations with integral operators 
by means of straightforward algebraic and differential operation on smooth 
functions. 



2 in plasma physics: an example 



In plasma physics pseudodifferential operators are found on dealing with 
linear wave propagation as a consequence of the dispersion combined with 
space- and time- variations of the background plasma. 

As an example of how symbol calculus can be applied, let us consider a 
cold isotropic plasma for which the constitutive relationship for electromag- 
netic waves is given in implicit form by the generalized Ohm's law which 
relates the induced current density J to the wave electric field E, namely, 

— J = -P- E 6 

dt 4vr y ' 

Wp e (r, t) being the squared electron plasma frequency which depends on 
both the point-location r and time t through the electron density profile. 
Equation © can be cast in the form Q J = E with Q the first order partial 
differential operator with total symbol q(uj,r,t) = = _ ^ 7T ( u '+ lu ) _ 

^pe ^pe 

= q + i + q , where v = v(r,t) / is a zero order symbol playing 

^pe 

the role of a fictitious collisional frequency. In the high frequency limit, 
the properties of the operator Q are dominated by the symbol q+\ which, 
however, is not unique being defined apart from the arbitrary function v (for 
sake of completeness let us mention that, for symbols that can be written 
as a sum of homogeneous functions, the order of homogeneity allows to 
distinguish uniquely the leading term [Tj; here, however, such a property is 
ignored as it does not hold in the general case). 
Let us search for a ^DO P such that 

PQ-I = R (7) 

where I is the identity operator (with symbol 1) and R a remainder of order 
— oo. A solution P of equation (J7J) is an approximate inverse (parametrix) of 
Q, thus, one gets J + R J = P E where the remainder R J can be neglected 
in the semiclassical limit with the result that J(r, t) ~ [P E)(r, t). 

The needed solution of equation is obtained on exploiting the 
product-rule ((Sa|) for left-reduced symbols which yields 



dp-i dq+i 

P-lQ+l +P-1Q0 +P-2Q+1 Qf ^ =1 > ( 8 ) 

where x = (r, ct) and k = (k, —u/c), = p~\ + p_2 + • • • being the left 
symbol of the parametrix P. Equation Q can be solved iteratively with the 



result that, to lowest significant orders, 
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(9a) 



(9b) 



Ait(uj + iv) 2 

In particular, one should note that p_i is the same as the conductivity of a 
stationary plasma with the time-dependence added adiabatically (adiabatic 
approximation) whereas the next order term p_2 explicitly accounts for the 
time variations of the plasma frequency. The induced current is related to 
the electric displacement by dt (D — E) = 4-7T J which can be solved for D on 
noting that the parametrix of dt corresponds to the symbol i(u> + iv)^ 1 — 
v(u + iv)~ 2 + • • • and applying again the product-rule (|5a|) with the result 
that D ~ £ E, the dielectric operator £ being a \&DO of order zero with 
symbol (cf., equations ©) 



e(u,r,t) = 1 



pe 



(u + ivY 
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It is worth noting that the symbol v(r, t) has no effects in the high frequency 
regime (as can be shown on Taylor-expanding with respect to v/uj <C 1), 
consistently with the considered generalized Ohm's law © which does not 
accounts for collisions; indeed, v has been introduced to deal the zero of q, 
though also a different (more formal) approach exists 1 which is based on 
smoothly cutting off the symbol near uj = 0. 

Finally, the resulting wave equation for the electric field reads 

^E)+Vx(VxE) = 0, (11) 

which follows directly from Maxwell's equations (with no external drive). 

The foregoing calculation of the response of a non-magnetized plasma 
exemplifies a general procedure that can be applied to more general cases 
(including, e.g., the effects of the magnetic field) and proves that the wave 
equation for electromagnetic waves in plasmas is, in general, a pseudodiffer- 
ential equation. 



3 Semiclassical asymptotics of ^DO 



A general solution of the inhomogeneous wave equation Q amounts to the 
sum of a solution of the corresponding homogeneous problem, D • -0 = 0, 
plus a particular solution of the whole inhomogeneous equation. 



The geometrical optics solution. Let us first address the homoge- 
neous problem and search for (asymptotic) solutions in the form (eikonal 
ansatz) [21 E] 

rj)(x) = e{x)A{x)e iS{ - x \ (12) 

e(x) and A{x) being the slowly varying (i.e., d x (e, A) ~ (e,A)/L with L 
defined in Q) polarization (unit) vector and complex amplitude, respec- 
tively, whereas the phase S(x), referred to as the real eikonal function, 
defines the local wavevector by k(x) = d x S(x) which is ordered according to 
\k(x)L\ ~ 5 _1 (^ +00 in the semiclassical limit). On substituting 1)12(1 into 
(^Q) with f = and expanding with respect to 5 one gets [2] 

,,, , , <N ,dd(k,x,x) d(eA) i d 2 S(x) d 2 d(k,x,x) , .» 
d(k, x, x) ■ (eA) - 1 — V • \ . ' r-^ — - ' • (eA) 

y ' v ; dh dx l 2dx l dx3 dhdkj v ; 

- i d2 f k %7 ] ■^ = ° {5 ' m+2 ^ (13) 

where k = k(x), i.e., the eikonal ansatz (|12j) allows to transform a pseudo- 
differential wave equation into a partial differential equation which has the 
form of an asymptotic series in 5 and thus can be solved iteratively. More 
specifically, with reference to a weakly non-Hermitian media for which d = 
dn+idA with dA = 0(5~ m+1 ), the lowest order approximation of (fl3|) leads 
to G(d x S(x), x) ■ e(x) = 0, D(k,x) = dn(k, x, x) being the local dispersion 
tensor. As far as mode conversion can be neglected jH], the polarization 
vector amounts to 

e(x) = e(d x S(x),x) (14a) 

with e(k, x) being an eigenvector of the dispersion tensor, corresponding to 
the eigenvalue D(k, x). Thereafter, the phase S(x) is determined by the first 
order partial differential equation (referred to as the eikonal equation) 

D(d x S(x),x) =0, (14b) 

which is analogous to the Hamilton-Jacobi equation for a Hamiltonian (me- 
chanical) system in the x-k phase space with D(k, x) as Hamiltonian function 
[HIS], so that the phase S(x) can be computed as the action relevant to such 
a Hamiltonian system. To the next order, equation I|1M)I yields the transport 
equation for the (squared) amplitude, namely, 



dx l 



— \A? 

dki 



2 7A |A| 2 , 7 A-ee* 



, , \(d' 1 d H d 2 d H 



2\dkidx l dkidx 



(14c) 



where k = k(x) and x' = x. Equation (|14c|) is a continuity equation for 
the flux vector ^|^4| 2 which points toward the velocity V % = ^ tangent 
to the trajectories of the equivalent mechanical system in the configuration 
(x) space, referred to as the geometrical optics (GO) rays, 7,4 being the 
absorption coefficient. As a consequence, equation (|14cj) can be reduced to 
an ordinary differential equations along the GO rays and solved in parallel 
to the Hamilton's equations with limited computational efforts. 

It should be noted that the form l)14j) of the GO equations holds inde- 
pendently on whether one makes use of the exact symbol d(k, x,x') or any of 
its reduced forms d^ q ' p \k, qx + px'); in the latter case the dispersion tensor 
amounts to D(fc, x) = d^' p \k, x), whereas the absorption coeffcient reads 

which, in particular, takes the simplest 



7 J' p) = ee 



d (9 ' p) - (n l ) d 2p 
U A \y 2>dk i dx i 



form if the Weyl symbol is used, namely, p = \ and 7^' 2 = e* • d^ 2 ' 2 • e. 
However, the polarization vector, the Hamiltonian function and the absorp- 
tion coefficient are different depending on the choice of the asymptotic form 
of the wave operator (J2J). One can prove that different formulations yield 
GO solutions which differ for 0(5)-corrections [I]. 

The Wigner-Weyl kinetic formalism. The GO solution discussed 
above exhibits two limitations: it has a local validity and does not apply to 
the inhomogeneous wave equation. To some extent, such limitations can be 
overcome on representing the wavefield in the x-k phase space. With this 
aim, the wave equation (^Q) is written in the quadratic form [2] 

D • (^) = (fft) • (Dt)-i, (15) 

where (D^) _1 can be obtained by means of symbol calculus whereas (tpip^) 
and (fT^) are projection operators. On assuming that both (i^ip^) and (fT^) 
are ^DO, one can apply the product-rule (|5b|) for the corresponding Weyl 
symbols with the result that equation (|15|) is transformed into jl] 

d m . s + i{ d w )S } -1 — = n • (dWt)-i + . . . ) (16) 

S and N being the Weyl symbols of (if)^) and (fft), respectively. Equation 
(|16|) can be solved iteratively as for the corresponding equation (fT3)l in the 
GO method. In particular, for weakly non-Hermitian (d^^ = 0(<5~ m+1 )) 
and weakly driven (N • (d^^) -1 = 0(5~ m+1 )) media and far from mode 



conversion regions, one has S(k,x) = W(k,x)e(k,x)e* (k, x) with 2 



D(k,x)W(k,x) = 0, (17a) 
[W ,D}-»tf>W + (17b) 

D(k,x) and = 7^ 5 ' 5 being, respectively, the eigenvalue of d^ 1 ^ cor- 
responding to the eigenvector e and the absorption coefficient in the Weyl 
form (cf., comments after equation (|14c[) l . Equation ()17b|) has the form of 
a kinetic equation in the x-k phase space for the scalar function W(k,x), 
referred to as the Wigner function, with the constraint (j!7a[) representing 
the local dispersion relation D(k,x) = 0. The Wigner function is related to 
the wavefield intensity through the identity |t/?(x)| 2 = J £\ W(k,x) [2], 
however, one can show that the approximations underlying the asymptotic 
series (|TK|) are such that the solution of the equations (|17|) allows to com- 
pute the wavefield intensity averaged over the large scalelength l{^> 2ir/\k\), 
(|?/>(x)| 2 )£, rather than the exact intensity \xj){x)\ 2 

Let us note that the Poisson brackets in l|17b[) . which are due to the 
specific product-rule (|5b|) for Weyl symbols, endows the x-k space with the 
same Hamiltonian structure as in the GO technique, provided that the Weyl 
formulation is used. 

Paraxial solutions of pseudodifferential wave equations. In the 
GO method, both the polarization vector e and the complex amplitude A 
have been assumed to be slowly varying on the scalelength L of the medium 
nonuniformity. Such an ansatz is justified for the polarization e far from 
mode conversion regions (cf., equation (|14ajl 1 . but not for the amplitude 
A that, being transported along the GO rays, can develope its own scale- 
length w, called beamwidth in analogy to the case of wavebeams, namely, 
d x A ~ A/w. Depending on the geometry of the bundle of GO rays, the 
beamwidth w can be shorter that L, e.g., near focal points, and should be ex- 
plicitly accounted for in the asymptotic analysis. With this aim, quasi-optics 
methods [£] have been developed, for the case of monochromatic electromag- 
netic wavebeams propagating in stationary, spatially nondispersive media for 
which the relevant wave equation exhibits local operators only. Among such 
quasi-optics techniques, the complex geometrical optics (CGO) method jH] 
is proved to be a benchmark on the basis of which other quasi-optics so- 
lutions can be derived [Zj; in particular, within the paraxial approximation 
(2ir/\k\) < w( ~ ^/2irL/\k\) < L, the CGO method reduces to the beam 
tracing (BT) technique (Sj which has the advantage of being particularly 
suited for numerical analysis. 



On combining the CGO method with the paraxial approximation, one 
can prove that the beam tracing method can be used to solve generic pseu- 
dodifferential wave equations provided that the range of nonlocal effects is 
finite Let us describe the general features of such BT solutions on the 
basis of the simple analytically-tractable case addressed in section 2; in par- 
ticular, let us assume that the plasma is homogeneous and stationary with 
e(u},r,t) = n 2 (u>), n(uj) being the refractive index, and search for solution 
of in the form E = y t/j(x,t), the relevant equation for the scalar field 
t/}(x,t) amounting to 

with N the 'I'DO corresponding to the refractive index n(u>). The wavefield 
at the initial time t = is assumed to be a wave train characterized by the 
wavelength 2ir/ko and by a Gaussian envelope with the initial e _1 -width wq, 

2/2 

namely, ip(x,t = 0) = tpo(x) = A§e~ x ' w ° cos(kox). It is worth noting that 
the exact solution of such a problem could be formally achieved by means of 
the Fourier transform, but the explicit calculation of the wavefield requires 
an integration in both the frequency lo and wavevector k which, in gen- 
eral, is not straightforward both analytically and numerically. On the other 
hand, the BT method relys on a set of ordinary differential equations [5-8] 
the coefficients of which are readily obtained as derivatives of the symbol 
D(k ,lu) — — ^n^iuf) + A; 2 corresponding to the ^DO in (|18[) , The solu- 
tions of such equations yield the wavefield directly in the (x, t) space-time, 
namely, [3] 



ip(x,t) = A 



Wq 



w(t) 

X COS 



S (t) + k (x- v g t) + ^y(^ - v g tf - <p(t) . (19) 



which is shown in figure 1. Let us note that the BT solution (|19j) is local- 
ized around the space-time curve (x = v g t, t) with v g = v v j (l + jj^t-j) being 
the group velocity (v p = c/n, the phase velocity), whereas the e _1 -width 
w(t) = \l + 4 (( x r )^g t ) 1 2 Wq increases throughout the propagation so that 
the Gaussian envelope spreads out depending on the adimensional paxame- 
ter r = J [l + ^ + (^) 2 + ^ 0] which controls the dispersion effects 

(for a non dispersive medium r = 1 and the BT solution (|19|) reduces to the 
well-known solution of the d'Alembert wave equation). As a consequence 



(a) 



(b) 



Figure 1: Time evolution of a Gaussian wave train according to the BT solu- 
tion l|19|) (r = 0.2 and v g /v p — 0.8): (a) two constant-time sections showing the 
propagation and spreading of the wavefield and (b) the plot of the solution in the 
two-dimensional space-time showing the localization of the wavefield around the 
curve (x — v g t, t). 



conservation of the squared phase-averaged amplitude integrated with re- 
spect to x. The short-scale oscillations of the wavefield are accounted for by 
the cos-term in (jl9|) . the argument of which exhibits a rather fine structure: 
So(t) = —icot', with t' = t — Vgt/vp the retarded time, describes the temporal 
oscillations, the term ko(x — v g t) is just the initial phase k$x shifted along the 
space-time curve (x = v g t,t) whereas the term 2 ( x ~ v gt) 2 accounts for 
the curvature of the phase-fronts (viewed as curves in the two-dimensional 

k 2 w A 

space-time) with R(t) = {\ — r)v g t+ i ^ 1 "_ r ^ v t the radius of curvature; finally, 
the phase-shift ip(t) = arctan [ 2 ^~a^ g^] is due to the effect of the envelope 
spreading (the analogous of the Gouy shift for wavebeams E] ) ■ In equa- 
tion ([TO]) any function of the frequency should be evaluated at u = ujq, with 
u>o solution of the dispersion equation D(ko, u>o) = 0. 

Let us note that several features of the BT solution (|19|) are general. In 
particular, the wavefield can be written in form of a com^lex-eikonal wavelet 
[6-8] t/j(x, t) = $te{A(x, t)e %s ( x ^ } where both the complex amplitude A(x, t) 
and the comp/ex-eikonal function S have been expanded in Taylor series with 




respect to (x — v g t) in a neighbourhood of the space-time curve (x = v g t,t), 
with time t labelling the position along the curve and (x — v g t) being its 
distance from the considered point (x,t) (paraxial expansion), namely, 



A(x,t) = A (t)+O(5*), A (t) = A c 



wo 



S(x, t) = So(t) + k (x - Vgt) + 



1 



+ 



w(t) 
2i 



2 e -Mt) 



R(t) w 2 (t) 



(x-v g tf + 0(55), 



which holds within the paraxial approximation with \x — v g t\ < w(t) 
As a consequence, the whole wavefield around the reference curve (which 
in the general case amounts to a GO ray referred to as the reference ray) 
is completely determined by a set of parameters, e.g., the amplitude Ao(t), 
the phase So(t), the width w[t) and the curvature radius R(t), that satisfy 
ordinary differential equations along the reference ray, so that one can find 
paraxial solutions of rather complicated wave equations (including nonlocal 
effects) by simple computational means. 



4 The intrinsic Hamiltonian structure 

It has been stressed that in the paraxial solution (|19|) the wavefield is local- 
ized around a geometrical optics (GO) reference ray, which represents the 
trajectory of the whole wave train in the space-time. On the other hand, the 
reference ray depends on the adopted GO formulation (cf., comments after 
equation ()14c|) ) . so that different forms of the dispersion tensor Y)(k,x) cor- 
responds to different reference rays which can deflect one from the other for 
long-enough propagation. Hence, one should wonder what GO formulation 
yields the reference ray around which the wavefield is actually localized (let 
us note that the symbol of the operator (|TH|) is uniquely defined, as follows 
from equation (JIJ, thus, in this case, there is no ambiguity in the definition 
of the reference ray) . 

With this aim, it is worth noting that the symbol d of the wave operator 
(J2J) is usually obtained by means of symbol calculus (e.g., as in section 2) 
and thus it amounts to an asymptotic series, namely, 

d(k, x, x') ~ d (k, x, x') + d 1 (k, x, x') H (20) 

with di being of order m — I. Moreover, the leading term in the reduced 
symbol ©, i.e., d { q,p \k,x) = d (k, x,x), is an Hermitian tensor for weakly 
non-Hermitian media and does not depend on the adopted (q,p)- formulation; 



for this reason, do is referred to as the principal symbol of the operator. 
Therefore, in spite of separating the Hermitian and anti-Hermitian parts of 
the total symbol d, one should make use of the asymptotic expansion ()2U|) 
into both equations (|13|) and with the result that the dispersion tensor 
is determined by the unique principal symbol, namely, D(fc, x) = do(k, x, x), 
and yields the same polarization vector and Hamiltonian function in any GO 
formulation as well as in the Wigner-Weyl formalism. As for the absorption 
coefficient, one gets (cf., equation (|14c|l l 



7A = ee 



\ ( dd d 2 d 

A + 



2\dkjdx l dkjdx 



ee 



1\ d 2 B 



2/ dkidx 1 



for any (q,p), where the asymptotic series (J3J has been accounted for 
toghether with (J2UJ). 

In conclusion, on separating the principal symbol it is possible to identify 
an intrinsic dispersion tensor which does not depend on the considered form 
of the operator and thus yields the same polarization vector, Hamiltonian 
function and absorption coefficient in every formulations of the geometrical 
optics method as well as in the Wigner-Weyl kinetic formalism 5 . However, 
such a separation is, in general, defined apart from the transformation D — > 
D + g with g a lower-order Hermitian symbol, e.g., the principal symbol 
q+i of the operator Q given in section 2 is defined on adding and arbitrary 
lower-order symbol v. Such an arbitrarity should be dealt with in term of a 
Gauge transformation 3 . 
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